suppressWarnings({
  suppressPackageStartupMessages({
    library(data.table)
    library(magrittr)
    library(data.table)
    library(dplyr)
    library(ggplot2)
    library(lubridate)
    library(vars)
    library(scales)
    library(forecast)
    library(tseries)
  })
})

options(repr.plot.width = 20, repr.plot.height = 20, repr.plot.res = 500)

# Duomenu uzkrovimas

url <- "http://web.vu.lt/mif/a.buteikis/wp-content/uploads/2023/09/Macro_2023_SVAR.R"
source(url, encoding = "UTF-8")
DT_0 <- get_data(2012527)
df <- as.data.frame(DT_0)

[1] "Country: Italy"
[1] "Surastas failas 'duomenys_atsiskaitymas_2.RDS'. Užkraunami duomenys iš direktorijos..."

I dalis: Duomenų paruošimas

1. Prieš pradedant analizuoti duomenis - detaliau susipažinkite su minėtu straipsniu ir aprašykite:

a. Į kokį klausimą bandoma atsakyti (t. y., kodėl atliekamas ekonometrinis modeliavimas)?

Skaičiuojamos fiskalinių daugiklių reikšmės tam tikriems makroekonominiams rodikliams tam, kad suprasti (išsiaiškinti) kaip šios reikšmės gali kisti priklausomai nuo valstybės skolos (angl. public debt), ekonomikos augimo greičio (angl. pace of economic growth) bei skirtumo tarp esamo bei potencialaus BVP (angl. output gap).

Ekonometrinis modeliavimas (konkrečiai, SVAR) yra atliekamas dėl to, jog fiskalinių rodiklių reikšmės yra kintančios. Pastarosios kinta dėl įvairių priežasčių kur viena iš jų galėtų būti egzogeniniai fiskaliniai šokai (kurie atsiranda dėl fiskalinės politikos pokyčių, ekonominių krizių ir kt.).

Taip pat, viena iš esminių SVAR savybių yra ta, jog šis modelis gali įvertinti kiekvieno iš kintamųjų atsaką tuo pačiu metu. SVAR gali padėti įvertinti sudėtingus ryšius tarp skirtingų kintamųjų tuo pačiu metu.

Ekonometrinio modeliavimo tikslas yra įvertinti fiskalinių (makroekomominių rodiklių) daugiklių reikšmes, atsirandančias dėl vyriausybės išlaidų bei mokestinių pajamų fiskalinio šoko.

b. Kokie straipsniuose naudojami duomenys VAR modeliavimui - kokie kintamieji, kokio dažnumo, koks laikotarpis, kokios šalies(-ių) duomenys, kokie analizuojamų rodiklių matavimo vienetai? Jeigu kai kurie rodikliai matuojami valiuta - ar tai realios, ar nominalios kainos?

Duomenys: Eurozonos šalių makroekonominiai rodikliai nuo EMU (European Monetary Union) sukūrimo. Laikotarpis nuo 2000 m. iki 2016 m. (ketvirtiniai duomenys).

Kintamieji:

GDP (liet. BVP)
Primary Expenditure (liet. valstybės išlaidos)
Income and Wealth Taxes (liet. valstybės pajamos iš mokesčių, taikomų pajamoms ir turtui)
Production and Imports Taxes (liet. valstybės pajamos iš mokesčių, taikomų produkcijai ir importui) \ Debt (liet. valstybės skola)

Visi kintamiej (išskyrus GDP) yra pateikiami realiomis kainomis (Eur), vienam gyventojui (angl. per capita). Taip pat, "Debt" yra % nuo "GDP".

c. Ar laiko eilučių duomenys kaip nors transformuojami (pvz. pašalinamas sezoniškumas, logaritmuojama, imami skirtumai, skaičiuojama vienam gyventojui tenkanti rodiklio dalis ar pan.) ar laiko eilutės stacionarios?

Taip, laiko eiličių duomenys yra transformuojami. Straipsnyje minima, jog visi rodikliai (apart GDP) yra paskaičiuoti vienam gyvenotojui bei logaritmuoti. Taip pat visi kintamieji (įskaitant GDP) yra paimti su skirtumais vienetinės šaknies testo atžvilgiu.

d. Kam naudojamas valstybės skolos nuo BVP dalies rodiklis?

Šis rodiklis yra naudojamas kaip pseudo kintamasis (angl. dummmy variable) dideliam valstybės skolos lygiui su 60 % slenksčiu (angl. threshold). Nustatytas slenkstis atskiria šalis, kur sumos yra žemiau arba virš šios reikšmės.

2. Iš turimų duomenų atsirinkite ir pasiruoškite reikalingus:

a. Išbrėžkite pradinių duomenų laiko eilučių grafikus. Kokias tendencijas pastebite? Ar pastebite kokių nors struktūrinių pasikeitimų?

DT_0 %>% filter(na_item != "Gross domestic product at market prices") %>%
  ggplot(aes(x = time, y = values)) +
  geom_line(aes(group = na_item, color = na_item), linewidth = 0.5) +
  facet_wrap(~ na_item, scales = "free", nrow = 5) +
  scale_y_continuous(labels = label_number(big.mark = ",")) +
  theme_bw() +
  theme(legend.position = "none")

[1] "General government: Current taxes on income, wealth, etc., receivable": valstybės pajamos iš mokesčių, taikomų pajamoms ir turtui (mln. Eur). Grafike gana aiškiai matomas vyraujantis sezoniškumas. Tai yra normalu, kadangi šio rodiklio duomenyse nėra pašalintas sezoniškumas. Taip pat, pastebima augimo tendencija, kadangi grafiko kreivė yra kylanti.

[2] "General government: Government consolidated gross debt": valstybės skola (% nuo BVP). Rodiklio duomenyse nėra pašalintas sezoniškumas, tačiau grafike to nematyti. Kaip ir prieš tai buvusiu atveju, matoma augimo tendencija. Tiesa ji yra kur kas aiškesnė, kadangi panašu, jog kreivė nėra įtakota sezoniškumo. Matomas pastovus skolos augimas. Maždaug nuo 2009 m. iki ~ 2012 m. skolos lygis kito itin mažai ir išliko lygiagretus X ašiai. Vėliau nuo maždaug 2013 m. matomas gana drastiškas augimas, kuris tęsėsi iki 2015 m. Nuo 2015 m. iki 2020 m. skolos lygis išliko maždaug vienodas, su nedideliais svyravimais aukštyn/žemyn. Pastebima, jog nuo 2020 m. prasidėjus COVID19 pandemijai, skolos lygis smarkiai ir drastiškai išaugo. Galiausiai, apie 2022 m. nukrito lygiai taip pat greitai kaip ir buvo pakilęs.

[3] "General government: Taxes on production and imports, receivable": valstybės pajamos iš mokesčių, taikomų produkcijai ir importui (mln. Eur). Kadangi rodiklio duomenyse nėra pašalintas sezoniškumas tai lygiai kaip ir pirmojo rodiklio atžvilgiu, itin matomas vyraujantis sezoniškumas. Grafike galima įžvelgti augimo tendenciją, kadangi grafiko kreivė yra kylanti.

[4] "General government: Total general government expenditure": valstybės išlaidos (mln. Eur). Rodiklio duomenyse nėra pašalintas sezoniškumas, tad grafike įžvelgiama sezoniškumo. Grafike taip pat matoma augimo tendencija ir ji buvo ryškiausia maždaug nuo 2020 m. Vėlgi, tai yra normalu, kadangi būtent tuo metu prasidėjo COVID19 pandemija ir valstybės išlaidos išaugo.

[5] "Total Population (Total), 15 years or over": populiacija >= 15 metų (tūkst.). Duomenyse nėra pašalintas sezoniškumas, tačiau kitu atveju, tai yra tiesiog populiacijos dydis šalyje. Pastebima stipri augimo tendencija nuo 2005 m. iki 2015 m. Nuo 2015 m. iki 2020 m. populiacijos lygis išliko maždaug vienodas. Nuo 2020 m. pastebima populiacijos mažėjimo tendencija, kuri yra gana drastiška.

# GDP grafikai

DT_0 %>% filter(na_item == "Gross domestic product at market prices") %>%
  ggplot(aes(x = time, y = values)) +
  geom_line(aes(color = interaction(s_adj, unit)), linewidth = 0.5) + labs(title = "Gross domestic product at market prices") +
  facet_wrap(s_adj ~ unit, scales = "free") +
  theme_bw() +
  theme(legend.position = "none")

Seasonally and calendar adjusted data:

Kadangi duomenyse pašalintas sezoniškumas, tai grafikai yra konkretesni. Labiau pastebimos tendencijos ir turimas aiškesnis vaizdas apie rodiklio kitimą.

[1] "Chain linked volumes (2015), million euro": BVP realiomis kainomis. Grafike sunku įžvelgti vien tik augimo ar kritimo tendenciją, tačiau tai vis tik labiau atspindi kritimą. Didžiausias šokas buvo patirtas 2020 m., kai rodiklio reikšmė nuo 425 tūkst. krito žemiau nei 375 tūkst.
[2] "Current prices, million euro": BVP nominaliomis kainomis. Bendrai grafike yra matoma augimo tendencija, tačiau kaip ir kitų rodiklių grafikuose, taip ir čia matomas gan stiprus kritimas, įvykęs 2020 m. Po jo vėl matoma augimo tendencija, kuri buvo gan sparti.
[3] "Price index (implicit deflator), 2015=100, euro": BVP defliatorius. Matoma vientisa augimo tendencija be didesnių sukrėtimų/šokų, kurie neigiamai paveiktų augimą.

Unadjusted data:

Šiuo atveju matome tuos pačius grafikus, kurie iš esmės atspindi tas pačias tendencijas kaip ir praeiti. Kadangi duomenyse nėra pašalintas sezoniškumas, tai jis grafikuose pastebimas gana stipriai.

b. Pasinaudokite [eurostat metodika] ir įsitikinkite, kad duomenyse santykis tarp nominalaus ir realaus BVP atitinka BVP defliatoriaus rodiklį.

### GDP (current prices & chain-linked volumes)

# Seasonally adjusted
    
GDP_nominal_adj <- sum(DT_0$values[DT_0$na_item == "Gross domestic product at market prices" &
                   DT_0$unit == "Chain linked volumes (2015), million euro" &
                   DT_0$s_adj == "Seasonally and calendar adjusted data"])

GDP_real_adj <- sum(DT_0$values[DT_0$na_item == "Gross domestic product at market prices" &
                                   DT_0$unit == "Current prices, million euro" &
                                   DT_0$s_adj == "Seasonally and calendar adjusted data"])

GDP_deflator_adj <- (GDP_nominal_adj / GDP_real_adj) * 100

GDP_deflator_reported_adj <- mean(DT_0$values[DT_0$na_item == "Gross domestic product at market prices" & 
                                                DT_0$unit == "Price index (implicit deflator), 2015=100, euro" & 
                                                DT_0$s_adj == "Seasonally and calendar adjusted data"], 
                                  na.rm = TRUE)

GDP_deflator_adj / GDP_deflator_reported_adj

# Duomenyse santykis tarp nominalaus ir realaus BVP beveik atitinka BVP defliatoriaus rodiklį (apskaičiuotas BVP defliatorius
# yra šiek tiek didesnis nei tas, kuris pateiktas duomenyse)

# Seasonally unadjusted
    
GDP_nominal_unadj <- sum(DT_0$values[DT_0$na_item == "Gross domestic product at market prices" &
                                      DT_0$unit == "Chain linked volumes (2015), million euro" &
                                      DT_0$s_adj == "Unadjusted data (i.e. neither seasonally adjusted nor calendar adjusted data)"])

GDP_real_unadj <- sum(DT_0$values[DT_0$na_item == "Gross domestic product at market prices" &
                                   DT_0$unit == "Current prices, million euro" &
                                   DT_0$s_adj == "Unadjusted data (i.e. neither seasonally adjusted nor calendar adjusted data)"])

GDP_deflator_unadj <- (GDP_nominal_unadj / GDP_real_unadj) * 100

GDP_deflator_reported_unadj <- mean(DT_0$values[DT_0$na_item == "Gross domestic product at market prices" & 
                                                DT_0$unit == "Price index (implicit deflator), 2015=100, euro" & 
                                                DT_0$s_adj == "Unadjusted data (i.e. neither seasonally adjusted nor calendar adjusted data)"], 
                                  na.rm = TRUE)

GDP_deflator_unadj / GDP_deflator_reported_unadj

# Duomenyse santykis tarp nominalaus ir realaus BVP beveik atitinka BVP defliatoriaus rodiklį (apskaičiuotas BVP defliatorius
# yra šiek tiek didesnis nei tas, kuris pateiktas duomenyse).

1.03984760701418

1.04023751823835

c. Apskaičiuokite reikalingų rodiklių reikšmes realiomis kainomis (žr. [Wiki 1], [Wiki 2] ir [“Inflation adjustment”]);

# Modeliavimui naudojami kintamieji: 

# "General government: Taxes on production and imports, receivable"
# "General government: Current taxes on income, wealth, etc., receivable"
# "General government: Total general government expenditure"
# "General government: Government consolidated gross debt"
# "Gross domestic product at market prices"

df <- subset(df, 
             na_item != "Total Population (Total), 15 years or over" & 
               !(na_item == "Gross domestic product at market prices" & 
                   s_adj == "Unadjusted data (i.e. neither seasonally adjusted nor calendar adjusted data)") & 
               unit != "Chain linked volumes (2015), million euro" & 
               unit != "Price index (implicit deflator), 2015=100, euro",
             select = -geo)

taxes_prod_imp <- data.frame(taxes_pr_imp = df$values[df$na_item == "General government: Taxes on production and imports, receivable"])
taxes_wealth <- data.frame(taxes_w = df$values[df$na_item == "General government: Current taxes on income, wealth, etc., receivable"])
gov_exp <- data.frame(exp = df$values[df$na_item == "General government: Total general government expenditure"])
gov_debt <- df$values[df$na_item == "General government: Government consolidated gross debt"] # % from GDP
gdp <- data.frame(gdp = df$values[df$na_item == "Gross domestic product at market prices"]) # current prices
gdp_def <- DT_0$values[DT_0$na_item == "Gross domestic product at market prices" & 
                         DT_0$unit == "Price index (implicit deflator), 2015=100, euro" &
                         DT_0$s_adj == "Seasonally and calendar adjusted data"]

taxes_prod_imp$time <- DT_0$time[DT_0$na_item == "General government: Taxes on production and imports, receivable"]
taxes_wealth$time <- DT_0$time[DT_0$na_item == "General government: Current taxes on income, wealth, etc., receivable"]
gov_exp$time <- DT_0$time[DT_0$na_item == "General government: Total general government expenditure"]



df2 <- data.frame(taxes_prod_imp = taxes_prod_imp, taxes_wealth = taxes_wealth, gov_exp = gov_exp,
                  gov_debt = gov_debt, gdp = gdp, gdp_def = gdp_def)

# Inflation adjustment
    
taxes_prod_imp$inf_adj <- (taxes_prod_imp$taxes_pr_imp / df2$gdp_def) * 100
taxes_wealth$inf_adj <- (taxes_wealth$taxes_w / df2$gdp_def) * 100
gov_exp$inf_adj <- (gov_exp$exp / df2$gdp_def) * 100
gdp$inf_adj <- (df2$gdp / df2$gdp_def) * 100

d. Jeigu reikia - pritaikykite kitas transformacijas, tokias kaip vienam gyventojui tenkanti rodiklio dalis, logaritmavimas, sezoniškumo šalinimas ir pan.

options(repr.plot.width = 15, repr.plot.height = 5, repr.plot.res = 500)

# Per capita adjustment
    
taxes_prod_imp$per_capita <- taxes_prod_imp$inf_adj / DT_0$values[DT_0$na_item == "Total Population (Total), 15 years or over"]
taxes_wealth$per_capita <- taxes_wealth$inf_adj / DT_0$values[DT_0$na_item == "Total Population (Total), 15 years or over"]
gov_exp$per_capita <- gov_exp$inf_adj / DT_0$values[DT_0$na_item == "Total Population (Total), 15 years or over"]
gdp$per_capita <- gdp$inf_adj / DT_0$values[DT_0$na_item == "Total Population (Total), 15 years or over"]

# Seasonality adjustment: 'taxes_prod_imp', 'taxes_wealth', 'gov_exp'

df_ts <- ts(taxes_prod_imp$per_capita, frequency = 4, start = c(year(min(taxes_prod_imp$time)), quarter(min(taxes_prod_imp$time))))
plot(df_ts, main = "Seasonally unadjusted TS, Taxes on Production & Import")
decomp <- stl(df_ts, s.window = "periodic")
taxes_prod_imp$sa_adj <- c(seasadj(decomp))

df_ts_adj_sa <- ts(taxes_prod_imp$sa_adj, frequency = 4, start = c(year(min(taxes_prod_imp$time)), quarter(min(taxes_prod_imp$time))))
plot(df_ts_adj_sa, main = "Seasonally adjusted TS, Taxes on Production & Import")

#####
    
df_ts <- ts(taxes_wealth$per_capita, frequency = 4, start = c(year(min(taxes_wealth$time)), quarter(min(taxes_wealth$time))))
plot(df_ts, main = "Seasonally unadjusted TS, Taxes on Income, Wealth, etc.")
decomp <- stl(df_ts, s.window = "periodic")
taxes_wealth$sa_adj <- c(seasadj(decomp))

df_ts_adj_sa <- ts(taxes_wealth$sa_adj, frequency = 4, start = c(year(min(taxes_wealth$time)), quarter(min(taxes_wealth$time))))
plot(df_ts_adj_sa, main = "Seasonally adjusted TS, Taxes on Income, Wealth, etc.")

#####

df_ts <- ts(gov_exp$per_capita, frequency = 4, start = c(year(min(gov_exp$time)), quarter(min(gov_exp$time))))
plot(df_ts, main = "Seasonally unadjusted TS, Government Expenditure")
decomp <- stl(df_ts, s.window = "periodic")
gov_exp$sa_adj <- c(seasadj(decomp))

df_ts_adj_sa <- ts(gov_exp$sa_adj, frequency = 4, start = c(year(min(gov_exp$time)), quarter(min(gov_exp$time))))
plot(df_ts_adj_sa, main = "Seasonally adjusted TS, Government Expenditure")

e. Ar nagrinėjami rodikliai turi vienetinę šaknį? Užrašykite nulinę hipotezę, aprašykite pasirinkto testo esmę (pvz. ADF testo) ir patikrinkite šią hipotezę kiekvienam rodikliui. Jeigu rodikliai turi vienetines šaknis - pritaikykite atitinkamas transformacijas į tai atsižvelgti.

# Testui naudojamas Augmented Dickey-Fuller (ADF) testas, kuris kiek pranašesnis nei įprastas DF testas tuo, jog naudoja 
# 'lagged' laiko eil. skirtumus. Tai padeda užtikrinti, jog atsitiktinės paklaidos nėra koreliuotos (angl. serially corellated). 
# Taip pat ADF testas yra lankstesnis, kadangi 'lagged' skirtumų įtraukimas leidžia turėti aukštesnės eilės serijinę koreliaciją.
# Taigi, ADF yra tinkamas įvairesniems duomenims nei, kad įprastas DF.
    
# H0: vienetinė šaknis egzistuoja
# H1: vienetinė šaknis neegzistuoja
    
adf.test(taxes_prod_imp$sa_adj, alternative = "stationary")
adf.test(taxes_wealth$sa_adj, alternative = "stationary")
adf.test(gov_exp$sa_adj, alternative = "stationary")
adf.test(gdp$per_capita, alternative = "stationary")

# taxes_prod_imp: p-value < 0.05 -> H0 atmetamas
# taxes_wealth: p-value > 0.05 -> H0 neatmetamas, vienetinė šaknis egzistuoja
# gov_exp: p-value > 0.05 -> H0 neatmetamas, vienetinė šaknis egzistuoja
# gdp: p-value > 0.05 -> H0 neatmetamas, vienetinė šaknis egzistuoja
    
# Kadangi kintamiesiems 'taxes_wealth', 'gov_exp' bei 'gdp' p-value buvo didesnė nei 0.05, tai reiškia, jog vienetinė šaknis
# egzistuoja. Kitaip tariant, šie laiko eil. duomenys nėra stacionarūs. Kintamiesiems bus atliekamos transformacijos ir tikrinama
# ar stacionarumas atsirado.
    
adf.test(diff(taxes_wealth$sa_adj), alternative = "stationary")

# Vieną kartą pritaikius 'lagged differences' funkciją rezultatas išliko nepakitęs ir duomenys toliau yra stacionarūs.
    
adf.test(diff(taxes_wealth$sa_adj, differences = 2), alternative = "stationary")

# Pritaikius diferencijavimą du kartus, atsirado stacionarumas. Teigiama, jog 'taxes_wealth' yra stac. laiko eil.
    
adf.test(diff(gov_exp$sa_adj), alternative = "stationary")

# Ta pati situacija kaip ir prieš tai, p-value > 0.05. Duomenys vis dar nestacionarūs.
    
adf.test(diff(gov_exp$sa_adj, differences = 2), alternative = "stationary")

# Atlikus papildomą transformaciją, p-value < 0.05. Teigiama, jog 'gov_exp' yra stac. laiko eil.
    
adf.test(diff(gdp$per_capita), alternative = "stationary")

# 'gdp' atveju užteko vieną kartą pritaikyti 'lagged differences' f-ją. Teigiama, jog 'gdp' yra stac. laiko eil.

Warning message in adf.test(taxes_prod_imp$sa_adj, alternative = "stationary"):
"p-value smaller than printed p-value"

	Augmented Dickey-Fuller Test

data:  taxes_prod_imp$sa_adj
Dickey-Fuller = -4.5454, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

	Augmented Dickey-Fuller Test

data:  taxes_wealth$sa_adj
Dickey-Fuller = -2.9726, Lag order = 4, p-value = 0.1793
alternative hypothesis: stationary

	Augmented Dickey-Fuller Test

data:  gov_exp$sa_adj
Dickey-Fuller = -3.4045, Lag order = 4, p-value = 0.06196
alternative hypothesis: stationary

	Augmented Dickey-Fuller Test

data:  gdp$per_capita
Dickey-Fuller = -3.3584, Lag order = 4, p-value = 0.06935
alternative hypothesis: stationary

	Augmented Dickey-Fuller Test

data:  diff(taxes_wealth$sa_adj)
Dickey-Fuller = -3.3567, Lag order = 4, p-value = 0.06975
alternative hypothesis: stationary

Warning message in adf.test(diff(taxes_wealth$sa_adj, differences = 2), alternative = "stationary"):
"p-value smaller than printed p-value"

	Augmented Dickey-Fuller Test

data:  diff(taxes_wealth$sa_adj, differences = 2)
Dickey-Fuller = -6.6211, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

	Augmented Dickey-Fuller Test

data:  diff(gov_exp$sa_adj)
Dickey-Fuller = -2.9225, Lag order = 4, p-value = 0.1999
alternative hypothesis: stationary

Warning message in adf.test(diff(gov_exp$sa_adj, differences = 2), alternative = "stationary"):
"p-value smaller than printed p-value"

	Augmented Dickey-Fuller Test

data:  diff(gov_exp$sa_adj, differences = 2)
Dickey-Fuller = -7.3776, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary

	Augmented Dickey-Fuller Test

data:  diff(gdp$per_capita)
Dickey-Fuller = -4.0099, Lag order = 4, p-value = 0.01431
alternative hypothesis: stationary

f. Galutinius atrinktus ir sutvarkytus rodiklius priskirkite atskiram kintamajam (kaip data.table, ar tibble tipo duomenų masyvą) ir išsibrėžkite tik galutinių atrinktų ir sutvarkytų rodiklių grafikus.

options(repr.plot.width = 20, repr.plot.height = 20, repr.plot.res = 500)

df_final <- data.frame(time = taxes_prod_imp$time, taxes_prod_imp = taxes_prod_imp$sa_adj, 
                       taxes_wealth = c(NA, NA, diff(taxes_wealth$sa_adj, differences = 2)),
                       gov_exp = c(NA, NA, diff(gov_exp$sa_adj, differences = 2)), 
                       gdp = c(NA, diff(gdp$per_capita)))

df_final_long <- melt(setDT(as.data.table(df_final)), measure.vars = 2:5, variable.name = "variable")

df_final_long %>%
  ggplot(aes(x = time, y = value)) +
  geom_line(aes(group = variable, color = variable), linewidth = 0.5) +
  facet_wrap(~ variable, scales = "free", nrow = 5) + # galima prideti nrow = 5; pabandyti pažiūrėti kitą kartą arba pridėti du grafikus: vienas su nrow ir kitas be
  theme_bw() +
  theme(legend.position = "none")

Warning message:
"Removed 5 rows containing missing values (`geom_line()`)."

Transformacijos, kurios buvo atliktos kintamiesiems:

• taxes_prod_imp: inflation adjustment, seasonality adjustment, per capita adjustment;
• taxes_wealth: inflation adjustment, seasonality adjustment, per capita adjustment, diferencijavimas (2x);
• gov_exp: inflation adjustment, seasonality adjustment, per capita adjustment, diferencijavimas (2x);
• gov_debt: -;
• gdp: inflation adjustment, per capita adjustment, diferencijavimas.

Labiausiai pasisekusi kintamojo transformacija buvo "taxes_prod_imp". Iš grafiko pastebima, jog sezoniškumo praktiškai nėra. Kas liečia kintamuosius "taxes_wealth" bei "gov_exp", sezoniškumas tam tikrose vietose dar vyrauja, tačiau grafikai atrodo kur kas geriau nei tie, kurie buvo išbrėžti naudojant pradinius (neapdorotus) duomenis.

Kas liečia "gov_debt", tai šiam kintamajam jokios transformacijos nebuvo atliekamos, kadangi duomenys yra procentiniai. Kintamajam "gdp" buvo atliktos tik dvi transformacijos, tačiau tiek pirminis tiek galutinis grafikas iš esmės atrodo identiškai. Pagrindinė to priežastis yra tai, jog pradiniuose duomenyse buvo pateikta "seasonally and calendar adjusted" GDP versija.

II dalis: Adekvataus VAR modelio sudarymas

3. Pasirinkite VAR modelio eilę, pagrįskite šį pasirinkimą. Sudarykite ir įvertinkite VAR modelį.

df_final <- na.omit(df_final)
rownames(df_final) <- NULL

VARselect(df_final[, 2:5], lag.max = 12, type = "const")$selection # 'const' - įtraukia 'intercept' komponentę, bet ne 'trend' ar 'seasonal'
VARselect(df_final[, 2:5], lag.max = 4, type = "const")$selection

# Pasirenkama, jog didžiausias lag skaičius būtų 12. Visos metrikos siūlo skirtingus lag kiekius. Lag skaičiaus pasirinkimas nėra
# pagrįstas. Pasirinkta atsitiktinai.
# Įdomumo dėlei papildomai pasirinktas straipsnyje naudojamas lag skaičius - 4.

AIC(n)

12

HQ(n)

3

SC(n)

2

FPE(n)

6

AIC(n)

4

HQ(n)

2

SC(n)

2

FPE(n)

4

VAR_mdl.1 <- VAR(df_final[, 2:5], p = 12, type = "const")
VAR_mdl.2 <- VAR(df_final[, 2:5], p = 6, type = "const")
VAR_mdl.3 <- VAR(df_final[, 2:5], p = 4, type = "const")
VAR_mdl.4 <- VAR(df_final[, 2:5], p = 3, type = "const")
VAR_mdl.5 <- VAR(df_final[, 2:5], p = 2, type = "const")

lapply(summary(VAR_mdl.1)$varresult, function(x){
  tibble::rownames_to_column(round(coef(x), 4) %>% as.data.frame(), "Coefficient")
})

$taxes_prod_imp

A data.frame: 49 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.4016	0.2379	1.6883	0.1256
taxes_wealth.l1	-0.1295	0.2747	-0.4715	0.6485
gov_exp.l1	-0.2698	0.1524	-1.7705	0.1104
gdp.l1	0.2497	0.1559	1.6016	0.1437
taxes_prod_imp.l2	0.2203	0.2697	0.8169	0.4351
taxes_wealth.l2	-0.4881	0.5805	-0.8409	0.4222
gov_exp.l2	-0.3123	0.2790	-1.1192	0.2920
gdp.l2	-0.0632	0.2061	-0.3065	0.7662
taxes_prod_imp.l3	-0.0036	0.3199	-0.0114	0.9912
taxes_wealth.l3	-1.0576	0.8406	-1.2581	0.2400
gov_exp.l3	-0.3583	0.3644	-0.9831	0.3512
gdp.l3	-0.0236	0.1375	-0.1719	0.8673
taxes_prod_imp.l4	0.1854	0.3344	0.5543	0.5929
taxes_wealth.l4	-1.0885	1.0156	-1.0718	0.3117
gov_exp.l4	-0.3445	0.3998	-0.8618	0.4112
gdp.l4	0.0213	0.0842	0.2527	0.8061
taxes_prod_imp.l5	-0.6420	0.3500	-1.8346	0.0998
taxes_wealth.l5	-0.8075	1.0015	-0.8063	0.4409
gov_exp.l5	-0.3730	0.4187	-0.8910	0.3961
gdp.l5	-0.0795	0.0687	-1.1575	0.2769
taxes_prod_imp.l6	-0.2767	0.3334	-0.8299	0.4280
taxes_wealth.l6	-0.8691	0.8811	-0.9865	0.3497
gov_exp.l6	-0.6920	0.4172	-1.6588	0.1315
gdp.l6	0.0832	0.0683	1.2180	0.2542
taxes_prod_imp.l7	0.1705	0.3513	0.4854	0.6390
taxes_wealth.l7	-0.3952	0.7629	-0.5180	0.6169
gov_exp.l7	-0.8650	0.4416	-1.9587	0.0818
gdp.l7	0.0228	0.0680	0.3352	0.7452
taxes_prod_imp.l8	0.6770	0.3593	1.8842	0.0922
taxes_wealth.l8	-0.0813	0.6991	-0.1163	0.9100
gov_exp.l8	-1.0150	0.4545	-2.2329	0.0524
gdp.l8	-0.0995	0.0767	-1.2970	0.2269
taxes_prod_imp.l9	-0.1680	0.3545	-0.4739	0.6469
taxes_wealth.l9	-0.0857	0.6631	-0.1293	0.9000
gov_exp.l9	-1.0565	0.4323	-2.4438	0.0371
gdp.l9	-0.0018	0.0811	-0.0227	0.9824
taxes_prod_imp.l10	0.1804	0.3248	0.5555	0.5921
taxes_wealth.l10	0.1425	0.5950	0.2395	0.8161
gov_exp.l10	-0.6749	0.3849	-1.7534	0.1134
gdp.l10	-0.0883	0.0799	-1.1048	0.2979
taxes_prod_imp.l11	-0.1837	0.3414	-0.5381	0.6036
taxes_wealth.l11	0.0081	0.4573	0.0177	0.9863
gov_exp.l11	-0.5999	0.3113	-1.9270	0.0861
gdp.l11	-0.0513	0.0699	-0.7342	0.4815
taxes_prod_imp.l12	-0.8611	0.3375	-2.5516	0.0311
taxes_wealth.l12	-0.1545	0.2361	-0.6547	0.5291
gov_exp.l12	-0.4252	0.2016	-2.1095	0.0641
gdp.l12	0.0404	0.0561	0.7204	0.4895
const	1.5422	0.5992	2.5737	0.0300

$taxes_wealth

A data.frame: 49 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.0114	0.2478	-0.0460	0.9643
taxes_wealth.l1	-1.6635	0.2862	-5.8128	0.0003
gov_exp.l1	-0.0809	0.1587	-0.5098	0.6224
gdp.l1	0.3405	0.1624	2.0967	0.0655
taxes_prod_imp.l2	-0.4684	0.2809	-1.6675	0.1298
taxes_wealth.l2	-2.0244	0.6047	-3.3479	0.0086
gov_exp.l2	0.2219	0.2906	0.7636	0.4646
gdp.l2	-0.5407	0.2147	-2.5180	0.0329
taxes_prod_imp.l3	-0.2927	0.3332	-0.8785	0.4025
taxes_wealth.l3	-2.3717	0.8756	-2.7085	0.0241
gov_exp.l3	0.0369	0.3796	0.0973	0.9246
gdp.l3	0.3007	0.1432	2.1003	0.0651
taxes_prod_imp.l4	-0.4664	0.3484	-1.3389	0.2134
taxes_wealth.l4	-1.8172	1.0579	-1.7178	0.1200
gov_exp.l4	0.0485	0.4164	0.1164	0.9099
gdp.l4	0.0233	0.0878	0.2659	0.7963
taxes_prod_imp.l5	0.2210	0.3645	0.6062	0.5594
taxes_wealth.l5	-1.2737	1.0432	-1.2209	0.2531
gov_exp.l5	-0.1388	0.4361	-0.3183	0.7575
gdp.l5	0.0056	0.0716	0.0776	0.9399
taxes_prod_imp.l6	0.7411	0.3473	2.1337	0.0616
taxes_wealth.l6	-1.1324	0.9178	-1.2339	0.2485
gov_exp.l6	-0.2592	0.4346	-0.5965	0.5655
gdp.l6	-0.0224	0.0712	-0.3154	0.7596
taxes_prod_imp.l7	0.1863	0.3659	0.5091	0.6229
taxes_wealth.l7	-0.6881	0.7947	-0.8659	0.4091
gov_exp.l7	-0.2805	0.4600	-0.6097	0.5572
gdp.l7	-0.0984	0.0708	-1.3900	0.1980
taxes_prod_imp.l8	0.1114	0.3743	0.2978	0.7726
taxes_wealth.l8	-0.9160	0.7282	-1.2579	0.2401
gov_exp.l8	-0.3140	0.4735	-0.6632	0.5238
gdp.l8	-0.0659	0.0799	-0.8249	0.4308
taxes_prod_imp.l9	-0.3863	0.3692	-1.0462	0.3227
taxes_wealth.l9	-0.9084	0.6908	-1.3151	0.2210
gov_exp.l9	-0.2096	0.4503	-0.4653	0.6527
gdp.l9	-0.0020	0.0845	-0.0234	0.9818
taxes_prod_imp.l10	-0.3198	0.3383	-0.9453	0.3692
taxes_wealth.l10	-0.6166	0.6198	-0.9949	0.3458
gov_exp.l10	-0.0496	0.4009	-0.1237	0.9043
gdp.l10	-0.0332	0.0832	-0.3983	0.6997
taxes_prod_imp.l11	-0.0597	0.3556	-0.1680	0.8703
taxes_wealth.l11	-0.5117	0.4764	-1.0740	0.3108
gov_exp.l11	-0.2183	0.3242	-0.6734	0.5176
gdp.l11	0.0108	0.0728	0.1485	0.8852
taxes_prod_imp.l12	-0.0145	0.3515	-0.0413	0.9680
taxes_wealth.l12	-0.2011	0.2459	-0.8179	0.4345
gov_exp.l12	-0.0330	0.2099	-0.1573	0.8785
gdp.l12	-0.0615	0.0585	-1.0526	0.3200
const	0.9017	0.6242	1.4446	0.1825

$gov_exp

A data.frame: 49 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.5699	0.5501	-1.0359	0.3273
taxes_wealth.l1	-0.1747	0.6353	-0.2750	0.7895
gov_exp.l1	-1.5430	0.3524	-4.3785	0.0018
gdp.l1	-0.2186	0.3606	-0.6063	0.5593
taxes_prod_imp.l2	0.7648	0.6236	1.2264	0.2512
taxes_wealth.l2	0.0339	1.3425	0.0253	0.9804
gov_exp.l2	-1.7315	0.6453	-2.6833	0.0251
gdp.l2	-0.1588	0.4767	-0.3330	0.7468
taxes_prod_imp.l3	0.2191	0.7397	0.2962	0.7738
taxes_wealth.l3	0.4886	1.9441	0.2513	0.8072
gov_exp.l3	-1.4213	0.8428	-1.6864	0.1260
gdp.l3	0.1206	0.3179	0.3795	0.7131
taxes_prod_imp.l4	-0.3569	0.7735	-0.4614	0.6554
taxes_wealth.l4	1.2805	2.3487	0.5452	0.5989
gov_exp.l4	-1.0548	0.9245	-1.1409	0.2834
gdp.l4	-0.0585	0.1948	-0.3002	0.7709
taxes_prod_imp.l5	-0.1122	0.8093	-0.1387	0.8928
taxes_wealth.l5	1.4791	2.3162	0.6386	0.5390
gov_exp.l5	-0.8218	0.9682	-0.8487	0.4180
gdp.l5	0.1677	0.1589	1.0553	0.3188
taxes_prod_imp.l6	-0.1483	0.7711	-0.1924	0.8517
taxes_wealth.l6	0.7668	2.0376	0.3763	0.7154
gov_exp.l6	-0.6960	0.9648	-0.7214	0.4890
gdp.l6	0.0779	0.1580	0.4933	0.6336
taxes_prod_imp.l7	0.2812	0.8124	0.3461	0.7372
taxes_wealth.l7	-0.0206	1.7643	-0.0117	0.9910
gov_exp.l7	-0.1802	1.0214	-0.1765	0.8638
gdp.l7	-0.0145	0.1572	-0.0920	0.9287
taxes_prod_imp.l8	-0.0857	0.8310	-0.1031	0.9202
taxes_wealth.l8	0.2577	1.6168	0.1594	0.8769
gov_exp.l8	0.1445	1.0512	0.1375	0.8937
gdp.l8	-0.0361	0.1775	-0.2036	0.8432
taxes_prod_imp.l9	-0.0609	0.8197	-0.0743	0.9424
taxes_wealth.l9	0.5994	1.5336	0.3909	0.7050
gov_exp.l9	0.4978	0.9998	0.4979	0.6305
gdp.l9	-0.0905	0.1876	-0.4823	0.6411
taxes_prod_imp.l10	-0.3029	0.7512	-0.4032	0.6962
taxes_wealth.l10	0.4865	1.3761	0.3535	0.7318
gov_exp.l10	0.1956	0.8902	0.2197	0.8310
gdp.l10	0.0144	0.1848	0.0776	0.9398
taxes_prod_imp.l11	0.0419	0.7895	0.0531	0.9588
taxes_wealth.l11	0.6160	1.0577	0.5824	0.5746
gov_exp.l11	0.1926	0.7199	0.2676	0.7951
gdp.l11	-0.0354	0.1616	-0.2194	0.8312
taxes_prod_imp.l12	0.2939	0.7805	0.3766	0.7152
taxes_wealth.l12	0.3440	0.5459	0.6301	0.5443
gov_exp.l12	0.0994	0.4661	0.2133	0.8358
gdp.l12	-0.0265	0.1298	-0.2039	0.8430
const	0.0480	1.3858	0.0346	0.9731

$gdp

A data.frame: 49 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.1140	0.4304	-0.2648	0.7971
taxes_wealth.l1	0.1488	0.4972	0.2993	0.7715
gov_exp.l1	-0.2888	0.2758	-1.0474	0.3222
gdp.l1	0.6565	0.2822	2.3268	0.0450
taxes_prod_imp.l2	-0.3828	0.4880	-0.7844	0.4529
taxes_wealth.l2	0.7994	1.0505	0.7610	0.4661
gov_exp.l2	-0.2817	0.5049	-0.5579	0.5905
gdp.l2	0.1077	0.3730	0.2886	0.7794
taxes_prod_imp.l3	0.5056	0.5788	0.8735	0.4051
taxes_wealth.l3	0.7441	1.5212	0.4892	0.6364
gov_exp.l3	-0.3795	0.6595	-0.5754	0.5791
gdp.l3	-0.1306	0.2487	-0.5250	0.6123
taxes_prod_imp.l4	0.3724	0.6052	0.6153	0.5536
taxes_wealth.l4	0.4488	1.8378	0.2442	0.8125
gov_exp.l4	-0.4666	0.7234	-0.6450	0.5350
gdp.l4	0.0118	0.1525	0.0775	0.9399
taxes_prod_imp.l5	-0.2480	0.6333	-0.3917	0.7044
taxes_wealth.l5	0.4101	1.8124	0.2263	0.8261
gov_exp.l5	-0.0144	0.7576	-0.0189	0.9853
gdp.l5	-0.0443	0.1244	-0.3566	0.7296
taxes_prod_imp.l6	-0.0716	0.6034	-0.1186	0.9082
taxes_wealth.l6	0.3073	1.5944	0.1927	0.8514
gov_exp.l6	0.1170	0.7549	0.1549	0.8803
gdp.l6	-0.0066	0.1236	-0.0535	0.9585
taxes_prod_imp.l7	-0.3821	0.6357	-0.6011	0.5626
taxes_wealth.l7	0.4496	1.3806	0.3257	0.7521
gov_exp.l7	0.1564	0.7992	0.1957	0.8492
gdp.l7	-0.0290	0.1230	-0.2355	0.8191
taxes_prod_imp.l8	-0.0491	0.6502	-0.0756	0.9414
taxes_wealth.l8	0.5290	1.2651	0.4181	0.6856
gov_exp.l8	-0.2458	0.8225	-0.2988	0.7719
gdp.l8	-0.0316	0.1389	-0.2276	0.8250
taxes_prod_imp.l9	0.1050	0.6414	0.1637	0.8736
taxes_wealth.l9	0.5617	1.2000	0.4681	0.6509
gov_exp.l9	-0.1040	0.7824	-0.1329	0.8972
gdp.l9	0.0706	0.1468	0.4809	0.6421
taxes_prod_imp.l10	0.3882	0.5878	0.6604	0.5255
taxes_wealth.l10	1.1154	1.0768	1.0359	0.3273
gov_exp.l10	0.2252	0.6965	0.3233	0.7539
gdp.l10	-0.0896	0.1446	-0.6196	0.5509
taxes_prod_imp.l11	0.2502	0.6178	0.4050	0.6950
taxes_wealth.l11	0.9704	0.8276	1.1726	0.2711
gov_exp.l11	-0.1663	0.5633	-0.2952	0.7745
gdp.l11	-0.0320	0.1264	-0.2528	0.8061
taxes_prod_imp.l12	-0.3102	0.6107	-0.5079	0.6237
taxes_wealth.l12	0.5005	0.4272	1.1717	0.2714
gov_exp.l12	-0.2667	0.3647	-0.7313	0.4832
gdp.l12	0.0512	0.1015	0.5037	0.6266
const	-0.0659	1.0844	-0.0608	0.9529

lapply(summary(VAR_mdl.2)$varresult, function(x){
  tibble::rownames_to_column(round(coef(x), 4) %>% as.data.frame(), "Coefficient")
})

$taxes_prod_imp

A data.frame: 25 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.4390	0.2006	2.1886	0.0347
taxes_wealth.l1	0.1021	0.1654	0.6175	0.5405
gov_exp.l1	-0.0093	0.1187	-0.0781	0.9382
gdp.l1	0.0705	0.0480	1.4689	0.1499
taxes_prod_imp.l2	0.2688	0.2243	1.1981	0.2381
taxes_wealth.l2	0.2760	0.2913	0.9476	0.3492
gov_exp.l2	0.0027	0.1958	0.0140	0.9889
gdp.l2	0.0271	0.0506	0.5355	0.5954
taxes_prod_imp.l3	0.0330	0.2236	0.1475	0.8835
taxes_wealth.l3	0.0045	0.3482	0.0128	0.9899
gov_exp.l3	-0.1333	0.2279	-0.5847	0.5621
gdp.l3	0.0255	0.0520	0.4898	0.6270
taxes_prod_imp.l4	0.3768	0.2184	1.7250	0.0924
taxes_wealth.l4	0.0597	0.3398	0.1757	0.8615
gov_exp.l4	-0.2545	0.2011	-1.2655	0.2132
gdp.l4	-0.0477	0.0525	-0.9087	0.3691
taxes_prod_imp.l5	-0.4897	0.2085	-2.3487	0.0240
taxes_wealth.l5	-0.0145	0.2877	-0.0503	0.9602
gov_exp.l5	-0.1995	0.1528	-1.3057	0.1993
gdp.l5	0.0151	0.0476	0.3170	0.7529
taxes_prod_imp.l6	-0.0450	0.2214	-0.2030	0.8402
taxes_wealth.l6	-0.2497	0.1600	-1.5607	0.1267
gov_exp.l6	-0.1406	0.1021	-1.3767	0.1765
gdp.l6	0.0636	0.0411	1.5470	0.1299
const	0.4933	0.1812	2.7222	0.0096

$taxes_wealth

A data.frame: 25 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.1503	0.1852	-0.8115	0.4220
taxes_wealth.l1	-1.5293	0.1527	-10.0134	0.0000
gov_exp.l1	0.0879	0.1096	0.8023	0.4273
gdp.l1	0.0044	0.0443	0.0991	0.9216
taxes_prod_imp.l2	-0.3097	0.2071	-1.4951	0.1429
taxes_wealth.l2	-1.5678	0.2689	-5.8295	0.0000
gov_exp.l2	0.2870	0.1808	1.5876	0.1204
gdp.l2	0.0507	0.0468	1.0845	0.2848
taxes_prod_imp.l3	-0.3294	0.2064	-1.5953	0.1187
taxes_wealth.l3	-1.6673	0.3215	-5.1858	0.0000
gov_exp.l3	0.1482	0.2104	0.7042	0.4855
gdp.l3	0.0906	0.0480	1.8883	0.0664
taxes_prod_imp.l4	0.0453	0.2017	0.2247	0.8234
taxes_wealth.l4	-1.3137	0.3138	-4.1867	0.0002
gov_exp.l4	0.0282	0.1857	0.1521	0.8799
gdp.l4	0.0388	0.0485	0.8004	0.4283
taxes_prod_imp.l5	0.1546	0.1925	0.8031	0.4268
taxes_wealth.l5	-0.9041	0.2657	-3.4030	0.0016
gov_exp.l5	-0.0485	0.1411	-0.3439	0.7328
gdp.l5	-0.0054	0.0439	-0.1226	0.9031
taxes_prod_imp.l6	0.3858	0.2044	1.8872	0.0666
taxes_wealth.l6	-0.5504	0.1477	-3.7259	0.0006
gov_exp.l6	0.0435	0.0943	0.4619	0.6467
gdp.l6	-0.0049	0.0380	-0.1283	0.8985
const	0.2406	0.1673	1.4381	0.1584

$gov_exp

A data.frame: 25 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.5052	0.2843	-1.7769	0.0834
taxes_wealth.l1	-0.1516	0.2345	-0.6465	0.5217
gov_exp.l1	-1.5544	0.1682	-9.2390	0.0000
gdp.l1	0.1567	0.0681	2.3027	0.0267
taxes_prod_imp.l2	0.4285	0.3180	1.3476	0.1856
taxes_wealth.l2	0.1053	0.4129	0.2551	0.8000
gov_exp.l2	-1.7428	0.2776	-6.2791	0.0000
gdp.l2	0.0790	0.0718	1.0999	0.2781
taxes_prod_imp.l3	0.5495	0.3169	1.7337	0.0909
taxes_wealth.l3	0.3292	0.4936	0.6671	0.5086
gov_exp.l3	-1.4272	0.3231	-4.4177	0.0001
gdp.l3	-0.0114	0.0737	-0.1548	0.8778
taxes_prod_imp.l4	0.1345	0.3096	0.4344	0.6664
taxes_wealth.l4	0.7635	0.4817	1.5851	0.1210
gov_exp.l4	-1.0232	0.2851	-3.5894	0.0009
gdp.l4	-0.0735	0.0745	-0.9871	0.3297
taxes_prod_imp.l5	-0.5626	0.2956	-1.9035	0.0644
taxes_wealth.l5	0.8968	0.4078	2.1990	0.0339
gov_exp.l5	-0.6144	0.2166	-2.8370	0.0072
gdp.l5	0.0288	0.0674	0.4269	0.6718
taxes_prod_imp.l6	-0.3949	0.3138	-1.2584	0.2157
taxes_wealth.l6	0.3056	0.2268	1.3477	0.1855
gov_exp.l6	-0.4251	0.1447	-2.9368	0.0055
gdp.l6	0.0917	0.0583	1.5732	0.1238
const	0.4142	0.2569	1.6126	0.1149

$gdp

A data.frame: 25 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.3612	0.7671	0.4708	0.6404
taxes_wealth.l1	0.2096	0.6325	0.3313	0.7422
gov_exp.l1	0.0174	0.4539	0.0382	0.9697
gdp.l1	-0.1575	0.1836	-0.8577	0.3963
taxes_prod_imp.l2	-0.4179	0.8578	-0.4872	0.6288
taxes_wealth.l2	0.9916	1.1139	0.8902	0.3788
gov_exp.l2	-0.1098	0.7488	-0.1466	0.8842
gdp.l2	-0.1587	0.1937	-0.8193	0.4176
taxes_prod_imp.l3	-1.0041	0.8550	-1.1744	0.2474
taxes_wealth.l3	0.6631	1.3315	0.4980	0.6213
gov_exp.l3	-0.2019	0.8716	-0.2317	0.8180
gdp.l3	0.1608	0.1987	0.8091	0.4234
taxes_prod_imp.l4	0.5533	0.8354	0.6624	0.5116
taxes_wealth.l4	0.4339	1.2995	0.3339	0.7403
gov_exp.l4	-1.1727	0.7690	-1.5249	0.1354
gdp.l4	-0.0155	0.2009	-0.0769	0.9391
taxes_prod_imp.l5	-0.8464	0.7974	-1.0615	0.2950
taxes_wealth.l5	-0.1253	1.1003	-0.1139	0.9099
gov_exp.l5	-0.8727	0.5843	-1.4936	0.1433
gdp.l5	0.1372	0.1819	0.7540	0.4554
taxes_prod_imp.l6	0.4604	0.8467	0.5438	0.5897
taxes_wealth.l6	-0.4689	0.6118	-0.7663	0.4481
gov_exp.l6	-0.7927	0.3905	-2.0300	0.0492
gdp.l6	0.1312	0.1573	0.8342	0.4092
const	1.0721	0.6930	1.5471	0.1299

lapply(summary(VAR_mdl.3)$varresult, function(x){
  tibble::rownames_to_column(round(coef(x), 4) %>% as.data.frame(), "Coefficient")
})

$taxes_prod_imp

A data.frame: 17 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.3378	0.1800	1.8766	0.0665
taxes_wealth.l1	-0.0232	0.1417	-0.1637	0.8707
gov_exp.l1	0.0264	0.0947	0.2793	0.7812
gdp.l1	0.0819	0.0458	1.7874	0.0801
taxes_prod_imp.l2	0.1360	0.2164	0.6286	0.5325
taxes_wealth.l2	0.0419	0.2092	0.2002	0.8421
gov_exp.l2	0.0411	0.1274	0.3228	0.7482
gdp.l2	0.0744	0.0492	1.5136	0.1365
taxes_prod_imp.l3	-0.0067	0.1920	-0.0347	0.9725
taxes_wealth.l3	-0.0966	0.2078	-0.4648	0.6441
gov_exp.l3	-0.0271	0.1154	-0.2348	0.8154
gdp.l3	0.0810	0.0473	1.7108	0.0934
taxes_prod_imp.l4	0.2296	0.1858	1.2357	0.2225
taxes_wealth.l4	0.0074	0.1359	0.0542	0.9570
gov_exp.l4	-0.0555	0.0758	-0.7325	0.4673
gdp.l4	0.0111	0.0379	0.2917	0.7717
const	0.3583	0.1595	2.2458	0.0293

$taxes_wealth

A data.frame: 17 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.3023	0.1802	-1.6772	0.0999
taxes_wealth.l1	-1.3483	0.1419	-9.5000	0.0000
gov_exp.l1	-0.0204	0.0948	-0.2153	0.8304
gdp.l1	0.0052	0.0459	0.1138	0.9099
taxes_prod_imp.l2	-0.1570	0.2166	-0.7248	0.4720
taxes_wealth.l2	-1.1706	0.2095	-5.5877	0.0000
gov_exp.l2	0.1688	0.1275	1.3234	0.1918
gdp.l2	0.0197	0.0492	0.4001	0.6908
taxes_prod_imp.l3	-0.1806	0.1923	-0.9390	0.3523
taxes_wealth.l3	-0.9090	0.2081	-4.3681	0.0001
gov_exp.l3	0.1695	0.1156	1.4660	0.1490
gdp.l3	0.0596	0.0474	1.2587	0.2141
taxes_prod_imp.l4	0.3829	0.1861	2.0576	0.0450
taxes_wealth.l4	-0.3870	0.1361	-2.8442	0.0065
gov_exp.l4	0.1233	0.0759	1.6237	0.1108
gdp.l4	-0.0610	0.0380	-1.6074	0.1144
const	0.3033	0.1598	1.8985	0.0635

$gov_exp

A data.frame: 17 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.6580	0.2725	-2.4145	0.0195
taxes_wealth.l1	-0.3166	0.2146	-1.4751	0.1466
gov_exp.l1	-1.1748	0.1434	-8.1949	0.0000
gdp.l1	0.2041	0.0693	2.9432	0.0050
taxes_prod_imp.l2	0.2239	0.3276	0.6834	0.4976
taxes_wealth.l2	-0.2534	0.3167	-0.7999	0.4276
gov_exp.l2	-1.0471	0.1928	-5.4300	0.0000
gdp.l2	0.1658	0.0744	2.2276	0.0305
taxes_prod_imp.l3	0.5346	0.2907	1.8388	0.0720
taxes_wealth.l3	-0.2860	0.3146	-0.9090	0.3678
gov_exp.l3	-0.6801	0.1748	-3.8914	0.0003
gdp.l3	0.0160	0.0716	0.2229	0.8246
taxes_prod_imp.l4	-0.2312	0.2814	-0.8217	0.4152
taxes_wealth.l4	-0.1340	0.2057	-0.6512	0.5180
gov_exp.l4	-0.3685	0.1148	-3.2098	0.0023
gdp.l4	0.0342	0.0574	0.5959	0.5540
const	0.1566	0.2415	0.6482	0.5199

$gdp

A data.frame: 17 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.1615	0.6567	0.2459	0.8068
taxes_wealth.l1	0.4326	0.5171	0.8365	0.4069
gov_exp.l1	-0.0753	0.3455	-0.2180	0.8283
gdp.l1	-0.1484	0.1671	-0.8882	0.3788
taxes_prod_imp.l2	-0.1858	0.7894	-0.2354	0.8149
taxes_wealth.l2	0.8857	0.7633	1.1605	0.2515
gov_exp.l2	-0.2177	0.4647	-0.4685	0.6415
gdp.l2	-0.0700	0.1794	-0.3903	0.6980
taxes_prod_imp.l3	-1.0516	0.7006	-1.5010	0.1398
taxes_wealth.l3	0.4674	0.7582	0.6165	0.5404
gov_exp.l3	-0.1319	0.4212	-0.3132	0.7554
gdp.l3	0.1585	0.1726	0.9179	0.3631
taxes_prod_imp.l4	0.1074	0.6780	0.1585	0.8748
taxes_wealth.l4	0.2826	0.4958	0.5701	0.5712
gov_exp.l4	-0.4478	0.2766	-1.6189	0.1119
gdp.l4	0.0255	0.1383	0.1847	0.8543
const	1.1535	0.5821	1.9817	0.0531

lapply(summary(VAR_mdl.4)$varresult, function(x){
  tibble::rownames_to_column(round(coef(x), 4) %>% as.data.frame(), "Coefficient")
})

$taxes_prod_imp

A data.frame: 13 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.4087	0.1692	2.4150	0.0191
taxes_wealth.l1	-0.0860	0.1065	-0.8072	0.4231
gov_exp.l1	0.0300	0.0788	0.3803	0.7052
gdp.l1	0.0609	0.0428	1.4225	0.1606
taxes_prod_imp.l2	0.0900	0.1855	0.4849	0.6297
taxes_wealth.l2	-0.0184	0.1447	-0.1273	0.8992
gov_exp.l2	0.0993	0.0904	1.0991	0.2766
gdp.l2	0.0619	0.0415	1.4899	0.1421
taxes_prod_imp.l3	0.1330	0.1626	0.8178	0.4170
taxes_wealth.l3	-0.1536	0.1004	-1.5300	0.1318
gov_exp.l3	0.0358	0.0693	0.5166	0.6076
gdp.l3	0.0451	0.0365	1.2356	0.2220
const	0.4349	0.1372	3.1694	0.0025

$taxes_wealth

A data.frame: 13 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.2294	0.1955	-1.1735	0.2457
taxes_wealth.l1	-1.1135	0.1231	-9.0483	0.0000
gov_exp.l1	0.0494	0.0910	0.5432	0.5893
gdp.l1	-0.0235	0.0495	-0.4759	0.6360
taxes_prod_imp.l2	-0.0188	0.2143	-0.0876	0.9305
taxes_wealth.l2	-0.7048	0.1671	-4.2172	0.0001
gov_exp.l2	0.2347	0.1044	2.2476	0.0287
gdp.l2	-0.0551	0.0480	-1.1476	0.2562
taxes_prod_imp.l3	-0.0340	0.1879	-0.1808	0.8572
taxes_wealth.l3	-0.4380	0.1160	-3.7767	0.0004
gov_exp.l3	0.1374	0.0801	1.7152	0.0920
gdp.l3	0.0222	0.0422	0.5262	0.6009
const	0.3343	0.1585	2.1089	0.0396

$gov_exp

A data.frame: 13 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.7686	0.2804	-2.7414	0.0083
taxes_wealth.l1	-0.3029	0.1765	-1.7164	0.0918
gov_exp.l1	-1.0389	0.1305	-7.9582	0.0000
gdp.l1	0.2102	0.0710	2.9619	0.0045
taxes_prod_imp.l2	0.0492	0.3074	0.1602	0.8733
taxes_wealth.l2	-0.2442	0.2397	-1.0188	0.3129
gov_exp.l2	-0.7128	0.1498	-4.7597	0.0000
gdp.l2	0.2164	0.0688	3.1436	0.0027
taxes_prod_imp.l3	0.6719	0.2695	2.4933	0.0158
taxes_wealth.l3	-0.1043	0.1663	-0.6273	0.5331
gov_exp.l3	-0.2752	0.1149	-2.3948	0.0201
gdp.l3	0.0162	0.0605	0.2669	0.7905
const	0.0576	0.2273	0.2534	0.8009

$gdp

A data.frame: 13 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.0938	0.6218	0.1508	0.8807
taxes_wealth.l1	0.1596	0.3914	0.4077	0.6851
gov_exp.l1	0.0278	0.2895	0.0960	0.9239
gdp.l1	-0.1625	0.1574	-1.0326	0.3064
taxes_prod_imp.l2	-0.3661	0.6817	-0.5371	0.5934
taxes_wealth.l2	0.4063	0.5315	0.7644	0.4479
gov_exp.l2	0.0906	0.3321	0.2727	0.7861
gdp.l2	-0.0255	0.1526	-0.1672	0.8678
taxes_prod_imp.l3	-0.8076	0.5976	-1.3513	0.1822
taxes_wealth.l3	0.0540	0.3688	0.1464	0.8841
gov_exp.l3	0.2456	0.2548	0.9640	0.3393
gdp.l3	0.1054	0.1342	0.7850	0.4359
const	1.2843	0.5041	2.5477	0.0137

lapply(summary(VAR_mdl.5)$varresult, function(x){
  tibble::rownames_to_column(round(coef(x), 4) %>% as.data.frame(), "Coefficient")
})

$taxes_prod_imp

A data.frame: 9 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.4778	0.1451	3.2929	0.0017
taxes_wealth.l1	-0.0284	0.0945	-0.3001	0.7651
gov_exp.l1	0.0554	0.0580	0.9543	0.3438
gdp.l1	0.0306	0.0383	0.7992	0.4274
taxes_prod_imp.l2	0.1775	0.1480	1.1990	0.2353
taxes_wealth.l2	0.1006	0.0866	1.1613	0.2502
gov_exp.l2	0.1221	0.0490	2.4899	0.0156
gdp.l2	0.0370	0.0354	1.0439	0.3008
const	0.4073	0.1282	3.1774	0.0024

$taxes_wealth

A data.frame: 9 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.1392	0.1976	-0.7043	0.4840
taxes_wealth.l1	-0.9306	0.1287	-7.2332	0.0000
gov_exp.l1	0.0305	0.0790	0.3866	0.7005
gdp.l1	-0.0877	0.0522	-1.6809	0.0981
taxes_prod_imp.l2	-0.0406	0.2016	-0.2012	0.8412
taxes_wealth.l2	-0.2690	0.1180	-2.2797	0.0263
gov_exp.l2	0.2222	0.0668	3.3269	0.0015
gdp.l2	-0.0603	0.0483	-1.2489	0.2166
const	0.2147	0.1746	1.2299	0.2236

$gov_exp

A data.frame: 9 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.7534	0.2595	-2.9031	0.0052
taxes_wealth.l1	-0.4043	0.1690	-2.3931	0.0199
gov_exp.l1	-0.8168	0.1038	-7.8709	0.0000
gdp.l1	0.1774	0.0685	2.5897	0.0121
taxes_prod_imp.l2	0.6285	0.2647	2.3738	0.0209
taxes_wealth.l2	-0.3444	0.1549	-2.2227	0.0301
gov_exp.l2	-0.5309	0.0877	-6.0527	0.0000
gdp.l2	0.1089	0.0634	1.7181	0.0910
const	0.1527	0.2292	0.6659	0.5081

$gdp

A data.frame: 9 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.2294	0.5211	0.4402	0.6614
taxes_wealth.l1	0.1601	0.3392	0.4720	0.6387
gov_exp.l1	-0.1955	0.2084	-0.9381	0.3520
gdp.l1	-0.1547	0.1376	-1.1246	0.2653
taxes_prod_imp.l2	-1.0334	0.5316	-1.9441	0.0567
taxes_wealth.l2	0.4501	0.3111	1.4467	0.1533
gov_exp.l2	-0.1037	0.1761	-0.5889	0.5582
gdp.l2	0.0698	0.1273	0.5484	0.5855
const	0.9581	0.4603	2.0814	0.0417

Išbandžius visų metrikų pasiūlytus 'lag' parametro kiekius, gaunami tokie rezultatai (stat. reikšm. kint. kiekio atžvilgiu): Lag = 12: itin daug kintamųjų ir praktiškai visi stat. nereikšm. Pasirinkimas atmetamas.

Lag = 6: taxes_prod_imp = 3, taxes_wealth = 6, gov_exp = 8, gdp = 1

Lag = 4: taxes_prod_imp = 1, taxes_wealth = 6, gov_exp = 7, gdp = 1

Lag = 3: taxes_prod_imp = 2, taxes_wealth = 5, gov_exp = 7, gdp = 1

Lag = 2: taxes_prod_imp = 3, taxes_wealth = 3, gov_exp = 8, gdp = 1

Kai lag sk. didžiausias, stat. nereikšm. kintamųjų gaunasi itin daug. Turint omenyje, jog viso turima 70 stebinių, tiek daug parametrų gali sukelti peroptimizavimą (overparametrized/overfitted). Su mažesniais 'lag' skaičiais kintamųjų parametrų kiekis yra daugiau mažiau panašus. Panašu, jog optimaliausias variantas būtų, kai lag = 3.

4. Atlikite Grangerio priežastingumo analizę (papildomas šaltinis su pavyzdžiu: [Ch. 11.4.1, Example 68]) ir kiekvienam rodikliui nustatykite, kurie rodikliai yra Grangerio priežastis. Uždėkite VAR modelio statistiškai nereikšmingiems koeficientams nulinius paribojimus. (Pastaba: tolimesnėse užduotyse nebūtina naudoti VAR modelio su tiesiniais apribojimais.)

restrictions <- matrix(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1,
                         0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 0, 1,
                         1, 0, 1, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0,
                         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1),
                       nrow = 4, ncol = 13, byrow = TRUE)

VAR_mdl_restricted.4 <- restrict(VAR_mdl.4, method = "manual", resmat = restrictions)
lapply(summary(VAR_mdl_restricted.4)$varresult, function(x){
  tibble::rownames_to_column(round(coef(x), 4) %>% as.data.frame(), "Coefficient")
})

$taxes_prod_imp

A data.frame: 2 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	0.5445	0.1040	5.2346	0
const	0.5394	0.1234	4.3717	0

$taxes_wealth

A data.frame: 5 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_wealth.l1	-1.0495	0.0964	-10.8821	0.0000
taxes_wealth.l2	-0.7226	0.1348	-5.3594	0.0000
gov_exp.l2	0.1374	0.0570	2.4111	0.0189
taxes_wealth.l3	-0.5054	0.1039	-4.8626	0.0000
const	-0.0002	0.0074	-0.0316	0.9749

$gov_exp

A data.frame: 7 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
taxes_prod_imp.l1	-0.8002	0.2025	-3.9526	0.0002
gov_exp.l1	-1.0406	0.0956	-10.8881	0.0000
gdp.l1	0.2051	0.0670	3.0623	0.0033
gov_exp.l2	-0.6993	0.1125	-6.2157	0.0000
gdp.l2	0.2583	0.0611	4.2249	0.0001
taxes_prod_imp.l3	0.8006	0.2020	3.9639	0.0002
gov_exp.l3	-0.3389	0.0930	-3.6453	0.0006

$gdp

A data.frame: 1 × 5

Coefficient	Estimate	Std. Error	t value	Pr(>|t|)
<chr>	<dbl>	<dbl>	<dbl>	<dbl>
const	0.0029	0.0233	0.1267	0.8995

H0: kintamasis nėra kitų kintamųjų Granger priežąstis
H1: kintamasis yra kitų kintamųjų Granger priežąstis

causality(VAR_mdl_restricted.4, cause = "taxes_prod_imp", vcov. = vcovHC(VAR_mdl_restricted.4))

# Kadangi p-value < 0.05, tai H0 atmetamas. Teigiama, jog 'taxes_prod_imp' yra kitų kintamųjų Granger priežąstis
    
causality(VAR_mdl_restricted.4, cause = "taxes_wealth", vcov. = vcovHC(VAR_mdl_restricted.4))

# Kadangi p-value < 0.05, tai H0 atmetamas. Teigiama, jog 'taxes_wealth' yra kitų kintamųjų Granger priežąstis

causality(VAR_mdl_restricted.4, cause = "gov_exp", vcov. = vcovHC(VAR_mdl_restricted.4))

# Kadangi p-value < 0.05, tai H0 atmetamas. Teigiama, jog 'gov_exp' yra kitų kintamųjų Granger priežąstis

causality(VAR_mdl_restricted.4, cause = "gdp", vcov. = vcovHC(VAR_mdl_restricted.4))

# Kadangi p-value < 0.05, tai H0 atmetamas. Teigiama, jog 'gdp' yra kitų kintamųjų Granger priežąstis

$Granger

	Granger causality H0: taxes_prod_imp do not Granger-cause taxes_wealth
	gov_exp gdp

data:  VAR object VAR_mdl_restricted.4
F-Test = 6.786, df1 = 9, df2 = 216, p-value = 1.429e-08


$Instant

	H0: No instantaneous causality between: taxes_prod_imp and taxes_wealth
	gov_exp gdp

data:  VAR object VAR_mdl_restricted.4
Chi-squared = 11.452, df = 3, p-value = 0.009519

$Granger

	Granger causality H0: taxes_wealth do not Granger-cause taxes_prod_imp
	gov_exp gdp

data:  VAR object VAR_mdl_restricted.4
F-Test = 7.014, df1 = 9, df2 = 216, p-value = 7.02e-09


$Instant

	H0: No instantaneous causality between: taxes_wealth and taxes_prod_imp
	gov_exp gdp

data:  VAR object VAR_mdl_restricted.4
Chi-squared = 6.1574, df = 3, p-value = 0.1042

$Granger

	Granger causality H0: gov_exp do not Granger-cause taxes_prod_imp
	taxes_wealth gdp

data:  VAR object VAR_mdl_restricted.4
F-Test = 6.8373, df1 = 9, df2 = 216, p-value = 1.218e-08


$Instant

	H0: No instantaneous causality between: gov_exp and taxes_prod_imp
	taxes_wealth gdp

data:  VAR object VAR_mdl_restricted.4
Chi-squared = 6.5618, df = 3, p-value = 0.08726

$Granger

	Granger causality H0: gdp do not Granger-cause taxes_prod_imp
	taxes_wealth gov_exp

data:  VAR object VAR_mdl_restricted.4
F-Test = 7.3511, df1 = 9, df2 = 216, p-value = 2.466e-09


$Instant

	H0: No instantaneous causality between: gdp and taxes_prod_imp
	taxes_wealth gov_exp

data:  VAR object VAR_mdl_restricted.4
Chi-squared = 9.7078, df = 3, p-value = 0.02122

5. Atlikite įvertinto VAR modelio liekanų analizę. Ką šios analizės rezultatai pasako apie modelio adekvatumą?

Autokoreliacijos testas (serial autocorrelation)

H0: liekanos nėra autokoreliuotos
H1: liekanos yra autokoreliuotos

serial.test(VAR_mdl_restricted.4, lags.pt = 4, type = "PT.asymptotic")

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl_restricted.4
Chi-squared = 100.65, df = 16, p-value = 2.609e-14

$serial

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl_restricted.4
Chi-squared = 100.65, df = 16, p-value = 2.609e-14

serial.test(VAR_mdl_restricted.4, lags.pt = 6, type = "PT.asymptotic")

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl_restricted.4
Chi-squared = 128.18, df = 48, p-value = 3.144e-09

$serial

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl_restricted.4
Chi-squared = 128.18, df = 48, p-value = 3.144e-09

serial.test(VAR_mdl_restricted.4, lags.pt = 12, type = "PT.asymptotic")

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl_restricted.4
Chi-squared = 176.78, df = 144, p-value = 0.03287

$serial

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl_restricted.4
Chi-squared = 176.78, df = 144, p-value = 0.03287

serial.test(VAR_mdl_restricted.4, lags.pt = 13, type = "PT.asymptotic")

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl_restricted.4
Chi-squared = 183.34, df = 160, p-value = 0.09978

$serial

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl_restricted.4
Chi-squared = 183.34, df = 160, p-value = 0.09978

Atlikus autokoreliacijos testą VAR modeliui su tiesiniais apribojimais, autokoreliacija egzistuoja iki pat 13-to lag'o. T. y., tik nuo 13-to lag'o p-value > 0.05. Kas indikuoja, jog autokoreliacijos nėra (H0 neatmetamas). Neatrodo, jog šį VAR modelį su tiesiniais apribojimais būtų adekvatu naudoti toliau.

serial.test(VAR_mdl.4, lags.pt = 5, type = "PT.asymptotic")

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 48.113, df = 32, p-value = 0.03358

$serial

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 48.113, df = 32, p-value = 0.03358

serial.test(VAR_mdl.4, lags.pt = 6, type = "PT.asymptotic")

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 59.847, df = 48, p-value = 0.1173

$serial

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 59.847, df = 48, p-value = 0.1173

serial.test(VAR_mdl.4, lags.pt = 13, type = "PT.asymptotic")

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 118.07, df = 160, p-value = 0.9946

$serial

	Portmanteau Test (asymptotic)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 118.07, df = 160, p-value = 0.9946

Atlikus autokoreliacijos testą tiesiškai neapribotam VAR modeliui, autokoreliacija liekanose dingsta ties 6-tu lag'u. Panašu, jog tai indikuoja labiau adekvatų VAR modelį turimiems duomenims. Toliau bus naudojamas VAR modelis be tiesinių apribojimų ('VAR_mdl.4').

options(repr.plot.width = 15, repr.plot.height = 8, repr.plot.res = 500)

mdl_serial_test <- serial.test(VAR_mdl.4, lags.pt = 6, type = "PT.asymptotic")
residual <- residuals(VAR_mdl.4)

plot(mdl_serial_test, name = "taxes_prod_imp")

hist(residual[, "taxes_prod_imp"], probability = TRUE, main = "Histogram and EDF") # papildomai nubrėžiama histograma, kadangi grafikų panelėje ji nelabai matosi
lines(density(residual[, "taxes_prod_imp"]), col = "blue")

plot(mdl_serial_test, name = "taxes_wealth")

hist(residual[, "taxes_wealth"], probability = TRUE, main = "Histogram and EDF")
lines(density(residual[, "taxes_prod_imp"]), col = "blue")

plot(mdl_serial_test, name = "gov_exp")

hist(residual[, "gov_exp"], probability = TRUE, main = "Histogram and EDF")
lines(density(residual[, "taxes_prod_imp"]), col = "blue")

plot(mdl_serial_test, name = "gdp")

hist(residual[, "gdp"], probability = TRUE, main = "Histogram and EDF")
lines(density(residual[, "taxes_prod_imp"]), col = "blue")

Pastebima, jog visų kintamųjų ACF bei PACF grafikuose esančios vertikalios linijos ties tam tikru kiekiu lag'ų neišlenda aukščiau ar žemiau pasikliovimo intervalo (mėlyna punktyrinė linija). Tai reiškia, jog nėra jokios žymios autokoreliacijos. Kitaip tariant modelio kintamųjų liekanos yra baltasis triukšmas (angl. white noise).

Ta pati tendencija liečia ir squared ACF/PACF grafikus. Esančios vertikalios linijos neišlinda už pasikliovimo intervalo ribų (tik vienetiniai atvejai kur viena ar dvi linijos yra šiek tiek virš numatytos ribos).

Taip pat panašu, jog kintamųjų "taxes_prod_imp", "taxes_wealth" bei "gov_exp" liekanų histogramos tenkina normalųjį skirstinį (kintamojo "gdp" atveju tas nelabai matosi). Tiesa, "taxes_prod_imp" histograma yra pasislinkusi į kairę (angl. left skewed), "taxes_wealth" pasislinkusi į dešinę (angl. right skewed), o "gov_exp" iš visų kintamųjų liekanų yra pati simetriškiausia ir labiausiai atitinkanti normalųjį pasiskirstymą. Panašu, jog įvertintas modelis yra daugiau mažiau adekvatus pagal kintamųjų liekanų histogramas.

Homoskedastiškumo testas

H0: liekanos yra homoskedastiškos
H1: liekanos nėra homoskedastiškos

print(arch.test(VAR_mdl.4, lags.multi = 6, multivariate.only = TRUE

	ARCH (multivariate)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 610, df = 600, p-value = 0.3797

$arch.mul

	ARCH (multivariate)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 610, df = 600, p-value = 0.3797

Kadangi p-value > 0.05, tai H0 neatmetamas. Teigiama, jog liekanos yra homoskedastiškos. Tai yra viena iš indikacijų, jog įvertintas modelis yra adekvatus.

Liekanų normalumo testas

H0: liekanos tenkina normalųjį skirstinį
H1: liekanos netenkina normalaus skirstinio

print(normality.test(VAR_mdl.4, multivariate.only = TRUE))

$JB

	JB-Test (multivariate)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 36.823, df = 8, p-value = 1.241e-05


$Skewness

	Skewness only (multivariate)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 5.9717, df = 4, p-value = 0.2013


$Kurtosis

	Kurtosis only (multivariate)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 30.851, df = 4, p-value = 3.284e-06

$jb.mul
$jb.mul$JB

	JB-Test (multivariate)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 36.823, df = 8, p-value = 1.241e-05


$jb.mul$Skewness

	Skewness only (multivariate)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 5.9717, df = 4, p-value = 0.2013


$jb.mul$Kurtosis

	Kurtosis only (multivariate)

data:  Residuals of VAR object VAR_mdl.4
Chi-squared = 30.851, df = 4, p-value = 3.284e-06

Nors iš ankščiau aptartų grafikų buvo galima įžiūrėti, jog liekanos tenkina normalumo sąlygą, tačiau pagal testo rezultatus matoma, jog prie 'JB-Test (multivariate)' esantis p-value < 0.05, H0 atmetamas. Teigiama, jog liekanos netenkina normalumo sąlygos. Bendrasis liekanų normalumo testas ('JB-Test') rodo nenormalumą. Pagrindinė to priežastis yra eksceso koeficientas (angl. kurtosis). Pastarojo p-value < 0.05. Tuo tarpu asimterijos (angl. skewness) p-value > 0.05. Tai reikštų, jog liekanos turi "sunkesnes" / "lengvesnes" uodegas (angl. tail) nei, kad normaliojo skirstinio atveju.

Liekanų struktūrinio lūžio testas

H0: liekanose yra struktūrinis lūžis
H1: liekanose nėra struktūrinio lūžio

mdl_cumsum <- stability(VAR_mdl.4, type = "OLS-CUSUM")
plot(mdl_cumsum)

Kadangi nei vieno iš kintamųjų grafiko kreivė neišlenda už raudonos tiesės ribų, tai H0 atmetamas. Teigiama, jog liekanose nėra struktūrinio lūžio (t. y., koeficientai yra stabilūs per numatytą duomenų laiko tarpą).

Modelio adekvatumo vertinimas

Liekanų autokoreliacijos testas: liekanos nėra autokoreliuotos
Homoskedastiškumo testas: liekanos yra homoskedastiškos
Normalumo testas: liekanos netenkina normalumo sąlygos
Liekanų struktūrinio lūžio testas: liekanose nėra struktūrinio lūžio

Kadangi įvertintas VAR modelis tenkino tris iš keturių testų, tai galima teigti, jog jis yra ganėtinai adekvatus. Didžiausias minusas yra, jog sudaryto modelio liekanos netenkina normalumo sąlygos. Kaip buvo pastebėta ankščiau, toks testo rezultatas buvo įtakotas eksceso koeficiento, kurio p-value < 0.05. Kita vertus, asimetrijos p-value > 0.05.

6. Atlikite įvertinto modelio IRF analizę, bei FEVD analizę - pasirinkite vieną rodiklį, iš kurio ateina impulsas (pvz. ortogonalus impulsas ateina iš iš valstybės išlaidų) ir pakomentuokite grafikų rezultatus.

# IRF (Impulso - Atsako) analizė (IRF: Impulse-Response Function)

set.seed(123)
for(i in c("taxes_prod_imp", "taxes_wealth", "gov_exp")){
  plot(irf(VAR_mdl.4, impulse = "gdp", response = i, n.ahead = 20, boot = TRUE))
}

'taxes_prod_imp'

Laikotarpio pradžioje matomas teigiamas 'gdp' šokas, kuris indikuoja, jog padidėjus 'gdp, 'taxes_prod_imp' (valstybės pajamos iš mokesčių, taikomų produkcija ir importui) taip pat padidėja. Atsakas į šoką yra didžiausias apie 4-ąjį ketvirtį, vėliau jis gan ženkliai mažėja. Bėgant laikui tampa neaišku kaip impulsas iš 'gdp' įtakoja 'taxes_prod_imp'.

'taxes_wealth'

Laikotarpio pradžioje 'gdp' sukeltas šokas panašu, jog neigiamai įtakoja 'taxes_wealth', apie 4-tą ketvirtį matomas gana nemažas 'taxes_wealth' padidėjimas. Vėliau jis šokinėja gana agresyviai kol galiausiai bėgant laikui nublanksta ir nėra aišku kaip 'gdp' šokas toliau įtakoja 'taxes_wealth' reikšmes.

'gov_exp'

Laikotarpio pradžioje matomas gana ryškus padidėjimas 'gov_exp' kas atrodo gana logiška. Didesnis 'gdp', didesnės valstybės išlaidos ('gov_exp'). Panašu, jog atsakas didžiausias pirmaisiais periodais, o vėliau kaip ir ankstesniuose grafikuose jis nublanksta.

# FEVD analizė (Prognozės paklaidų dispersijos dekompozicija)
    
VAR_mdl.4_decomp <- fevd(VAR_mdl.4, n.ahead = 20)
plot(VAR_mdl.4_decomp)

taxes_prod_imp dispersija daugiausiai susideda iš nukrypimo nuo taxes_prod_imp. Kiti makrorodikliai panašu, jog įtakos daro itin mažai. Didžiausią įtaką daro gdp, kurio dalis vizualiai atrodo didžiausia.

taxes_wealth dispersija daugiausiai susideda iš nukrypimo nuo taxes_wealth, tačiau salyginai nemažą dalį (lyginant su likusiais rodikliais) sudaro nukrypimas nuo taxes_wealth_prod. Kas galėtų indikuoti, jog šie mokesčiai tarpusavyje susiję.

gov_exp dispersija daugiausiai susideda iš nukrypimo nuo gov_exp. Tuo tarpu nukrypimai nuo taxes_wealth ir taxes_prod_imp turi panašią įtaką. Tai yra logiška, kadangi pgr. valstybės pajamų šaltinis yra mokesčiai. Iš to seka kaip daug valstybė gali išleisti.

gdp dispersija daugiausiai susideda iš nukrypimo nuo gdp bei taxes_prod_imp. Kadangi antrasis rodiklis yra būtent taxes_prod_imp, tai galėtų indikuoti, jog šioje vietoje gaunami mokesčiai sudaro kur kas didesnę dalį bendrų pajamų, gaunamų iš mokesčių apskritai (lyginant taxes_prod_imp ir taxes_wealth).

III dalis: SVAR modelio specifikacija ir rezultatų analizė

7. Sudarykite SVAR modelį:

a. Koks SVAR modelis buvo sudarytas straipsnyje - A-modelis, B-modelis, ar AB-modelis?

Straipsnyje buvo sudarytas A-modelis. Diagonalės elementai yra vienetai, kas priklauso sudarinėjant standartinį SVAR modelį. Po diagonale esantys elementai žymi vienalaikius (angl. contemporaneous) vieno kintamojo efektus kitam kintamajam. Tuo tarpu nuliai matricoje žymi tokius kintamusuosius, kurie pastarojo efekto neturi.

b. Kiek nežinomųjų galime turėti, kad sistema SVAR modelyje būtų tiksliai identifikuota?

Viso nežinomųjų: $\frac{N(N + 1)}{2} = \frac{4(4 + 1)}{2} = \frac{4 \cdot 5}{2} = 10$, kur $N - kint. sk.$

Taigi, viso galima turėti 10 kintamųjų norint, jog sistema SVAR modelyje būtų tiksliai identifikuota.

c. Kokios straipsnyje įvardytos prielaidos, pagal kurias uždedami apribojimai?

Naudojama Cholesky dekompozicija (angl. Cholesky decomposition) ir kintamieji yra išdėstyti taip, jog pirmasis kintamasis ("Taxes on Income and wealth") neturi jokio vienalaikio efekto likusiems kintamųjų šokams (kurie yra žymimi "u"). Tokio tipo išdėstymą matricoje indikuoja pirmoje eilutėje esantys nuliai (išskyrus pirmąjį skaičių, kuris yra diagonalėje bei nekreipiant dėmesio į paskutinį skaičių, kuris yra įvertintas parametras). Tęsiant toliau, antrasis kintamasis ("Taxes on Production and imports") reaguoja tik į pirmąjį šoką ir nereaguoja į likusius (tą indikuoja nulis, esanatis 2-oje eil. 3-iame stulp.). Trečiasis kintamasis ("Primary expenditure") reaguoja į pirmus du šokus, tačiau ne į ketvirtąjį (tą indikuoja 3-ioje eil. 4-ame stulp. esantis nulis). Galiausiai, ketvirtasis kintamasis ("GDP") gali reaguoti į visus prieš tai buvusius šokus (tą indikuoja nenuliniai koeficientai 4-oje eil.).

d. Užrašykite šiuos apribojimus matriciniu pavidalu atitinkamoje SVAR modelio specifikacijoje ir įvertinkite šį modelį.

restrictions.SVAR <- matrix(c(1, 0, 0, 0,
                              NA, 1, 0, 0,
                              NA, NA, 1, 0,
                              NA, NA, NA, 1),
                            nrow = 4, ncol = 4, byrow = TRUE)

SVAR.mdl <- SVAR(VAR_mdl.4, Amat = restrictions.SVAR, max.iter = 10000)

print(SVAR.mdl)

SVAR Estimation Results:
======================== 


Estimated A matrix:
               taxes_prod_imp taxes_wealth gov_exp gdp
taxes_prod_imp         1.0000       0.0000  0.0000   0
taxes_wealth          -0.2974       1.0000  0.0000   0
gov_exp               -0.4211       0.5224  1.0000   0
gdp                   -2.1732       0.3375  0.1624   1

SVAR.mdl$LR$data.name
SVAR.mdl$LR$method
SVAR.mdl$LR$p.value
SVAR.mdl$LR$statistic

df_final[, 2:5]

'LR overidentification'

Chi^2: 0

Chi^2: 1374.99407124309

SVAR.mdl$A # įvertintų koef. matrica

A matrix: 4 × 4 of type dbl

	taxes_prod_imp	taxes_wealth	gov_exp	gdp
taxes_prod_imp	1.0000000	0.0000000	0.0000000	0
taxes_wealth	-0.2974250	1.0000000	0.0000000	0
gov_exp	-0.4211336	0.5224167	1.0000000	0
gdp	-2.1732305	0.3375159	0.1624419	1

SVAR.mdl$Ase # A paklaidų matrica

A matrix: 4 × 4 of type dbl

	taxes_prod_imp	taxes_wealth	gov_exp	gdp
taxes_prod_imp	0.0000000	0.0000000	0.0000000	0
taxes_wealth	0.1221694	0.0000000	0.0000000	0
gov_exp	0.1274586	0.1221694	0.0000000	0
gdp	0.1374509	0.1378361	0.1221694	0

8. Ištirkite įvertintą SVAR modelį:

a. Kaip atrodo šio modelio impulso-atsako funkcijų grafikai, naudojant impulsą iš to pačio rodiklio kaip 8 užduotyje? Ar jie skiriasi nuo VAR modelio IRF grafikų?

for(i in c("taxes_prod_imp", "taxes_wealth", "gov_exp")){
  plot(irf(SVAR.mdl, impulse = "gdp", response = i, n.ahead = 20, boot = TRUE))
}

par(mfrow = c(2, 2), mar = c(2.2, 2.2, 1.5, 0.5))
for(i in c("taxes_prod_imp", "taxes_wealth", "gov_exp")){
  V_T <- irf(VAR_mdl.4, impulse = "gdp", response = i, n.ahead = 20, boot = TRUE)
  S_T <- irf(SVAR.mdl, impulse = "gdp", response = i, n.ahead = 20, boot = TRUE)
  V_T <- unlist(V_T$irf, use.names = FALSE)
  S_T <- unlist(S_T$irf, use.names = FALSE)
  #
  plot.ts(V_T, ylim = c(min(V_T, S_T), max(V_T, S_T)), 
          main = paste0("Ortogonalus impulsas iš GDP, atsakas ", i))
  lines(ts(S_T), col = "blue")
  abline(h = 0, lty = 2, col = "red")
  legend("right", legend = c("VAR", "SVAR"), lty = 1, col = c("black", "blue"), cex = 0.8)
}

Warning message in SVAR(x = varboot, Amat = restrictions.SVAR, max.iter = 10000):
"Convergence not achieved after 10000 iterations. Convergence value: 1.23539440677334e-07 ."

Lyginant SVAR ir VAR modelių IRF grafikus matomas itin didelis skirtumas. Naudojant tą patį impulsą, likusių makroekonominių rodiklių atsakai yra ženkliai didesni (tiek teigiami tiek neigiami). T. y., pasiekiamos aukštesnės reikšmės Y ašyje, tačiau pačios kreivės tendencija išlieka identiška: pikai ir žemumos yra identiškose vietose ir taip pat vienodai nublanksta bėgant laikui).

b. Remdamiesi straipsnyje nurodyta metodika, apskaičiuokite fiskalinius daugiklius. Kaip jūsų konkrečios šalies modelio fiskaliniai daugikliai skiriasi nuo straipnio 2 lentelės rezultatų?

irf_taxes_prod_imp <- irf(VAR_mdl.4, impulse = "taxes_prod_imp", response = "gdp", n.ahead = 4, boot = TRUE) # 'n.ahead = 4' - kadangi duomenys ketvirtiniai, tai 'irf' f-ja vertins šoką metams į priekį
irf_taxes_wealth <- irf(VAR_mdl.4, impulse = "taxes_wealth", response = "gdp", n.ahead = 4, boot = TRUE)
irf_gov_exp <- irf(VAR_mdl.4, impulse = "gov_exp", response = "gdp", n.ahead = 4, boot = TRUE) 

taxes_prod_imp_multiplier <- sum(irf_taxes_prod_imp$irf$taxes_prod_imp[ , "gdp"])
taxes_wealth_multiplier <- sum(irf_taxes_wealth$irf$taxes_wealth[ , "gdp"])
gov_exp_multiplier <- sum(irf_gov_exp$irf$gov_exp[ , "gdp"])

fiscal_multipliers <- c("taxes_prod_imp" = taxes_prod_imp_multiplier, "taxes_wealth" = taxes_wealth_multiplier,
                        "gov_exp" = gov_exp_multiplier)
fiscal_multipliers

taxes_prod_imp

-0.0333203373390602

taxes_wealth

0.00748890525247068

gov_exp

-0.0082628195028223

Lyginant įvertintus rodiklius su tais, kurie pateikti straipsnyje matomas itin didelis skirtumas. Šiuo atveju įvertinti rodikliai yra tik vienos šalies - Italijos. Tuo tarpu straipsnyje autoriai naudoja Eurozonos šalių makroekonominius rodiklius nuo EMU (European Monetary Union) sukūrimo. Turint omenyje, jog duomenys sudaryti iš keliolikos šalių, tai galbūt todėl daugikliai taip žymiai skiriasi. Kiekvienos iš šalių makroekonominiai rodikliai gali skirtis itin drastiškai, kas gali įtakoti didesnius fiskalinius daugiklius.